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Regional  seismograms  from  the  Soviet  JVE  explosion  were  modeled.  The  explosion  part 
of  the  source  contained  considerable  RDP  overshoot  with  possible  spall  contribution.  The 
tectonic  release  part  was  composed  of  stress  relaxation  and  secondary  high  frequency  sources. 

Pn,  Sn,  Rayleigh  and  Love  waves  were  modeled. 

SsPmp  phases  were  discovered  in  regional  waves  and  modeled  using  ray  theory  and 
synthetic  seismogram  computations.  SsPmp  is  an  important  wave  in  broadband  data  for  event 
under  the  crustal  waveguide  but  also  occurs  for  crustal  events  under  low  velocity  surface  layers. 

When  structure  is  known,  SsPmp  can  be  used  as  an  added  constraint  on  location  of  regional 
events. 

Teleseismic  source  inversion  using  broadband  SV  waves  was  investigated.  It  was  found 
that  SV/SH  inversion  can  constrain  source  radiation  patterns  as  well  as  P/SH  inversion  although 
information  on  receiver  structures  is  very  important.  Smaller  earthquake  events  become 
accessible  for  waveform  inversion  using  SV  and  SH  waveforms. 

Theory  and  computational  methods  were  explored  for  the  problems  of  plane  wave 
propagation  in  plane  layered,  general  anisotropic  media  and  point  sources  in  heterogeneous 
media.  A  teleseismic  body  wave  inversion  package  was  constructed  for  routine  use. 
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Executive  Summary 

At  the  beginning  of  the  conception  of  this  work  in  1990,  the  direction  of 
research  in  the  seismic  verification  field  was  quickly  being  refocussed  onto  the 
issues  of  a  proposed  comprehensive  test  ban  (CTB).  Whereas  past  seismic 
research  had  been  concerned  with  yield  determinations  from  relatively  large 
underground  nuclear  explosions,  the  proposed  CTB  significantly  broadened  the 
issues  of  event  detection,  location  and  identification.  It  became  clear  that 
possible  violations  of  a  CTB  treaty  would  probably  be  associated  with  small 
events  which  could  only  be  detected  at  regional  distances.  The  related  problems 
of  event  detection,  location,  and  identification  literally  became  entwined  on  the 
skein  of  complex,  high  frequency  regional  wave  propagation.  Although 
verification  Seismology  had  made  tremendous  progress  in  teleseismic  wave 
propagation  for  detecting,  locating  and  characterizing  seismic  sources  of  all 
types,  the  interest  in  the  regional  wave  propagation  regime  represented  a 
significant  challenge  for  solving  the  usual  problems  for  smaller  events  at  regional 
distances. 

Work  on  this  grant  has  concentrated  on  the  problems  of  regional  wave 
propagation  although  attention  to  teleseismic  problems  was  also  given.  A  major 
effort  was  made  to  understand  the  excitation  and  propagation  of  regional  body 
and  surface  waves  from  the  1988  Joint  Verification  Experiment  explosion  in  the 
Former  Soviet  Union.  Broadband  data  written  by  the  explosion  on  instruments 
deployed  by  the  Natural  Resources  Defense  Council  were  (and  still  are)  unique 
in  displaying  spectacular,  clear  regional  waveforms.  Study  of  these  data  yielded 
provocative  clues  to  the  nature  of  both  the  seismic  source  and  how  waves 
propagated  in  the  near-regional  distance  range.  These  results  are  summarized 
below  and  were  published  by  Langston  (1995). 
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An  interesting  regional  phase  for  events  occurring  under  the  crustal 
waveguide  was  discovered  by  accident  while  reviewing  data  recorded  by  the 
IRIS  GSN.  It  turns  out  that  the  phase  SsPmp  can  be  the  largest  arrival  on  a 
regional  seismogram  and  complicates  the  interpretation  of  S  waves.  This  phase 
is  described  by  an  SV  wave  which  leaves  upwards  from  the  source,  converts  to  a 
P  wave  at  the  free  surface  and  then  reflects  from  the  Moho.  Its  large  size  is  due 
to  the  combination  of  a  large  conversion  coefficient  at  the  free  surface  and  to  the 
fact  that  the  resulting  P  wave  reflection  at  the  Moho  is  usually  post-critical.  An 
application  of  using  this  phase  for  earthquake  location  was  described  in  the  final 
report  of  a  previous  contract  (Ecker,  1992).  General  characteristics  and 
observations  of  SsPmp  are  described  by  Langston  (1996)  and  are  again 
summarized  below. 

Much  of  the  seismicity  in  Eurasia  is  intraplate  in  type  and  often  consists  of 
small  to  moderate  earthquakes  which  can  only  be  recorded  at  a  few  teleseismic 
stations.  It  is  often  these  very  earthquakes  which  may  be  of  the  greatest  interest 
in  regional  studies  because  they  may  be  plentiful  and  yield  high  signal-to-noise 
signals  at  regional  distances.  However,  the  paucity  of  teleseismic  data  makes  it 
difficult  to  obtain  good  source  parameters  such  as  source  orientation,  depth  and 
time  function.  These  are  parameters  which  are  very  useful  in  understanding 
regional  waves  since  the  structure  effect  and  the  source  effect  must  be  separated. 
A  study  by  Wang  and  Langston  (1995)  (see  below)  addressed  the  source 
inversion  problem  for  combined  inversion  of  SV  and  SH  waves  from  earthquake 
sources.  Because  shear  waves  are  more  efficiently  generated  from  earthquakes 
than  are  P  waves,  events  with  magnitudes  less  than  5.5  or  so  often  write 
seismograms  with  well-defined  S  waveforms  but  not  P  waveforms.  The 
governing  equations  for  the  response  of  a  dislocation  source  show  that  the 
horizontal  radiation  pattern  for  SV  and  P  waves  is  the  same.  Thus,  it  should  be 
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possible  to  use  SV  and  SH  waveforms  for  determining  source  parameters  in  the 
same  way  that  P  and  SH  waveforms  are  used.  However,  most  seismologists 
avoid  using  the  SV  waveform  because  of  severe  wave  propagation  effects  due  to 
the  crustal  waveguide  which  produces  Shear-coupled  PI  phases.  Wang  and 
Langston  (1995)  show  that  it  is  indeed  possible  to  use  broadband  SV  waves  in  a 
source  inversion  if  reasonable  receiver  structures  are  available. 

Two  other  topics  in  teleseismic  wave  propagation  were  addressed  in  an 
exploration  of  future  wave  propagation  studies.  In  the  first,  Daejin  Kang 
developed  the  equations  needed  to  use  the  reciprocal  theorem  for  constructing 
the  far-field  solution  for  a  point  source  embedded  in  complex,  3  dimensional 
elastic  media.  The  solution  was  checked  using  plane  layer  Thomson-Haskell 
theory.  The  reciprocal  theorem  makes  it  possible  to  perform  far-field 
computations  of  point  sources  in  heterogeneous  media  by  using  relatively  simple 
numerical  plane  wave  models.  The  direction  of  this  research  is  pointed  to  using 
the  theory  with  plane  wave  3D  finite  difference  calculations  for  study  of  the 
effect  of  lateral  heterogeneity  on  far  field  body  waves  from  earthquakes  and 
explosions.  The  second  theoretical  topic  considered  was  that  of  synthetic  plane 
wave  seismograms  for  layered  anisotropic  media.  Daejin  Kang  extended  the 
theory  of  Keith  and  Crampin  (1977a;b;c)  to  include  a  free  surface  above  a  stack  of 
layers  over  a  halfspace.  Synthetic  seismograms  for  incident  P,  SV  or  SH  waves 
can  be  computed  for  arbitrary  elastic  anisotropy.  This  technique  will  find  use  in 
interpretation  of  split  teleseismic  S  waves,  such  as  SKS,  and  modeling  of 
anisotropic  receiver  responses. 

Finally,  the  author  spent  time  writing  a  comprehensive  synthetic 
seismogram  and  inversion  software  package  for  analysis  of  teleseismic  body 
waves.  Certainly,  this  is  not  a  new  approach  since  the  author  has  considered 
these  problems  since  1974.  However,  the  package  incorporates  aspects  of  the 
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Joint  Seismic  Program  Center  (JSPC)  database  software  package  with  a  thorough 
synthetic  seismogram  calculator  which  makes  broadband  inversions  fast  and 
routine.  It  is  also  straightforward  to  use  the  package  for  regional  waveforms 
with  the  addition  of  a  suitable  wavenumber  integration  code.  The  system  was 
presented  at  an  IRIS  shortcourse  on  Moment  Tensor  Inversion  given  on  the 
Berkeley,  California,  campus  December  15  and  16, 1995.  The  tutorial  is  presented 
in  the  Appendix  and  the  software  will  be  publicly  available  from  the  IRIS  DMS 
March  1996. 
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SUMMARY  of 

Anatomy  of  Regional  Phases  and  Source  Qiaracterization  of  the  Soviet  Joint 
Verification  Experiment,  Underground  Nuclear  Explosion 

by  C.  A.  Langston 

Broadband  data  from -the  NRDC  deployment  of  stations  in  Kazakhstan  at 
Karkaralinsk  (KKL)  and  Bayanual  (BAY)  are  subjected  to  detailed  waveform 
modeling  to  understand  the  nature  of  regional  wave  propagation  and  to 
determine  the  nature  of  the  explosion  time  function  and  sources  of  low-  and 
high-frequency  "tectonic  release".  Minor  modification  of  a  crustal  model  inferred 
from  DSS  studies  in  the  region  is  very  effective  in  explaining  the  arrival  times, 
amplitudes,  and  waveshapes  of  the  Pn  and  PmP  arrivals  at  both  stations.  The 
SH-wave  data  are  remarkable,  however,  for  showing  a  large,  high-frequency 
arrival  after  the  expected  and  observed  arrival  of  SmS  from  the  explosion  source. 
These  secondary  S  arrivals  imply  a  secondary  source  1  to  a  few  kilometers  to  the 
south  of  the  shot  point.  Shear  waves  from  this  source  dominate  the  seismic  data 
at  high  frequency  but  are  minor  at  low  frequency  where  shear  waves  and  surface 
waves  from  the  inferred  relaxation  of  the  explosion  cavity  dominate  the  wave 
field.  The  secondary  S-wave  source  is  inferred  to  represent  explosion-driven 
faulting  with  no  net  moment.  Significant  overshoot  in  the  RDP  is  required  to 
match  the  P-wave  data.  It  is  found  that  von  Seggern-Blandford  and  Haskell 
source  functions  are  incapable  of  producing  the  observed  overshoot.  However, 
the  Mueller-Murphy  source  model  is  very  effective  in  matching  the  clear 
overshoot  effects  seen  in  the  P  waves.  Spall  effects  are  also  considered  but 
appear  to  be  minor  and  serve  to  require  even  more  overshoot  in  the  RDP.  The 
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excitation  and  dispersion  of  Rayleigh  waves  in  the  period  band  of  1  to  10  sec  is 
seen  to  be  a  sensitive,  but  nonunique,  function  of  both  Vs  and  Vp  in  the  upper 
crustal  layer.  High  Poisson  ratios  in  this  layer  are  needed  to  explain  the  Rayleigh 
waveform. 
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SUMMARY  of 

The  SsPmp  Phase  in  Regional  Wave  Propagation 

by 

Charles  A.  Langston 

Regional  broadband  waveforms  of  deep-  and  intermediate-depth  earthquakes 
written  at  College,  Alaska,  and  Matsushiro,  Japan,  show  large  SsPmp  phases  that 
may  be  the  largest  phase  seen  on  the  seismogram.  SsPmp  is  created  when  the 
direct,  upgoing  SV-wave  converts  to  a  P-wave  at  the  free  surface  and 
subsequently  becomes  trapped  in  the  crust  because  of  postcritical  P-wave 
reflection  at  the  Moho.  SsPmp  may  arrive  between  the  P-  and  S-wave  phases  or 
it  may  arrive  after  S,  depending  on  source  depth  and  distance.  Thus,  it  can  be 
used  as  an  additional  constraint  on  source  depth  and  location.  It  is  important  to 
recognize  the  existence  of  this  phase  since  it  is  easily  confused  wiht  the  S-wave 
arrival,  resulting  in  erroneous  S-wave  arrival  time  picks,  or  it  may  be  interpreted 
as  a  split  shear  wave,  incorrectly  implpng  shear-wave  anisotropy  in  the  medium. 
Such  a  wave  propagation  effect  will  also  be  important  for  a  crustal  source  for 
structures  that  contain  low-velocity  layers  near  the  surface. 
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SUMMARY  of 

Use  of  Broadband  Teleseismic  SV  Waves  in  Earthquake  Source  Studies 

by 

Mingguang  Wang  and  Charles  A.  Langston 

A  generalized  inverse  technique  utilizing  a  moment  tensor  formalism  is  applied 
to  teleseismic,  broadband  SV,  SH,  and  P  wave  data.  The  SV  wave  Green's 
functions  are  generated  using  the  SKBJ  spectral  method,  which  allows  the  source 
and  the  receiver  regions  to  have  different  structures  and  accurately  simulates  the 
wave  propagation  complexities  of  SPL  and  S-coupled  P  waves  in  the  source  and 
receiver  areas.  P  and  SH  wave  Green’s  functions  are  calculated  by  the 
generalized  ray  method.  Comparison  of  joint  P/SH,  SV/SH  and  P/SV/SH 
source  parameter  inversions  of  the  February  21, 1991,  Bering  Sea  event  and  the 
April  25,  1992,  Cape  Mendocino,  California  event  demonstrate  that  the  SV 
waveform  can  effectively  take  the  place  of  the  P  waveform  in  teleseismic  source 
studies.  However,  receiver  and  source  structural  complexities  produce  poorer 
SV  waveform  fits,  compared  to  P  and  SH  fits.  Joint  inversion  of  all  waveform 
data  does  best,  however. 
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Teleseismic  Body  Wave  Radiation  from  a  Point  Source  in  Three-Dimensional 
Heterogeneous  Media  Using  the  Reciprocity  Theorem 

by 

Daejin  Kang 
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1.  Introduction 


It  is  generally  believed  that  scattered  energy  resulting  from  heterogeneity  near  the  source,  along  the 
raypath,  and  near  the  receiver  is  related  to  waveform  complexity.  Recently,  several  studies  have  been 
conducted  that  address,  in  detail,  the  effect  of  seismic  wave  propagation  through  heterogeneous  media 
(e.g.,  Vidale  and  Helmberger,  1988).  These  studies  found  that  structural  variations  along  the  propagation 
path  can  have  a  profound  effect  on  the  amplitude  and  duration  of  the  recorded  seismic  signal. 

We  focus  our  attention  on  the  source  region  for  these  types  of  study.  A  point  source  is  located  at 
either  in  or  near  three-dimensional  heterogeneous  medium.  Our  objective  is  the  determination  of  the  far- 
field  displacement  produced  by  this  source-medium  interaction  using  the  principle  of  seismic  reciprocity. 

The  reciprocity  theorem  is  a  useful  tool  to  study  teleseismic  radiation  from  seismic  sources  in  the 
Earth’s  crust  (White,  1960;  Gupta,  1965,  1966,  1967;  Ward,  1973;  Bouchon  ,  1976;  Mclaughlin  et  al., 
1987,  1988).  Though  most  of  the  results  could  have  been  obtained  through  solutions  to  the  forward  prob¬ 
lem  the  reciprocal  approach  to  the  problem  has  the  advantage  of  simplicity  and  may  provide  more  physical 
insight.  It  is  also  applicable  to  the  case  where  the  source  is  located  in  general  media  (heterogeneous, 
anisotropic).  In  the  forward  problem  the  complications  of  the  point  source  and  the  general  media  occur 
together,  making  the  solution  difficult.  On  the  other  hand,  the  approximate  solution  of  the  reciprocal 
problem  is  conceptually  far  easier  to  obtain.  In  the  reciprocal  problem  the  complications  of  the  point 
source  and  the  general  media  are  separated.  In  this  problem  the  point  source  occurs  in  a  homogeneous 
region  far  from  the  general  media.  The  wavefronts  are  approximately  plane  when  they  strike  the  hetero¬ 
geneous  and  anisotropic  zone.  This  response  of  the  complicated  media  is  more  easily  determined  for 
incident  plane  waves.  The  complexity  of  scattering  by  the  heterogeneous  and  anisotropic  zone  in  the 
near-field  of  a  point  source  is  avoided  by  considering  the  reciprocal  problem. 

In  the  present  study,  we  will  show  that  the  reciprocal  approach  to  the  problem  -  the  inversion  of 
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source  and  receiver  through  the  reciprocity  theorem  -  can  lead  to  the  same  results  in  a  much  simplified 
manner.  Examples  presented  in  this  proposal,  however,  are  restricted  to  the  flat  layers  with  explosion 
source. 

2.  Purposes 

The  purposes  of  the  proposed  research  are  ( 1)  to  learn  the  details  of  earthquake  and  explosion  sources 
(2)  to  determine  earth  structure  and  (3)  from  the  knowledge  of  sources  and  structures,  to  predict  the 
displacements  at  the  surface  of  the  earth  from  earthquakes  and  explosions.  In  particular,  we  will  focus  our 
attention  on  the  effects  of  laterally  varying  complex  structures  in  source  region  on  radiated  waveform  and 
spectrum. 

3.  Reciprocal  problem 

Reciprocity  is  expressed  as  a  symmetry  of  Green’s  tensor  for  linear  elastodynanucs.  This  result 
follows  from  the  fact  that  the  equations  of  motion  in  a  linear  heterogeneous,  anisotropic  elastic  medium 
with  homogeneous  boundary  conditions  form  a  self-adjoint  system.  The  theorem  has  been  shown  to  hold 
for  elastic  waves  (Knopoff  &  Gangi,  1959;  Gangi,  1970).  According  to  White  (1960),  the  theorem  of 
reciprocity  may  be  stated  as  follows  : 

If  in  a  bounded,  heterogeneous,  anisotropic,  elastic  medium,  a  transient  force  f(t)  applied  in  some  partic¬ 
ular  direction  a  at  some  point  P  (Fig.  1)  creates  at  a  second  point  Q  a  transient  displacement  whose 
component  in  some  direction  p  is  u(t),  then  the  application  of  the  same  force  f(t)  at  point  Q  in  direction 
p  will  cause  a  displacement  at  point  P  whose  component  in  the  direction  a  is  u(t). 

We  can  express  this  reciprocity  theorem  in  terms  of  Green’s  tensor.  Consider  the  Green  s  tensor 
Gni(x,t;^,0)  which  represents  the  n-th  component  of  displacement  at  position  x  and  time  t  caused  by  the 
application  of  an  instantaneous  force  of  unit  impulse  in  the  i  direction  at  position  ^  and  time  0.  The  reci¬ 
procity  principle  states  that  the  following  relation  holds  for  the  Green’s  tensor  (Aki  &  Richards,  1980) 
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G„i(x,t;^,0)  =  Gi„(^,t;x,0) 


(1) 


In  this  proposal,  we  wish  to  use  this  reciprocity  principle  to  simplify  the  determination  of  the  far-field 
displacement  from  an  explosive  point  source.  A  spherical  explosive  source  of  any  extent  is  equivalent  to 
three  equal  and  mutually  orthogonal  force  dipoles.  Therefore,  the  moment  tensor  representation  for  the 
explosive  source  is  given  by  (Geller,  1976) 

Un(x,t)  =  -  [Gnl,l  +  Gn2.2  +  GnS.s]  *  Mq  (2) 

Applying  the  reciprocity  theorem,  equation  (1),  to  the  Green’s  tensor  in  equation  (2),  we  find  that 

Un(x,t)  =  -  [Gin,l  +  G2n.2+  G3n.3]  *  Mq  0) 

This  equation  represents  the  n-th  component  of  displacement  at  Q  in  the  forward  problem  (Fig.  2)  is 
obtained  from  the  divergence  of  displacements  at  P  in  the  reciprocal  problem  in  which  the  source  is  a 
single  point  force  at  Q  directed  along  the  n-th  axis. 

4.  Teleseismic  body  wave  radiation  from  an  explosive  source  in  flat  layered  media. 

Now,  a  simple  example  is  presented  which  combines  the  reciprocity  theorem  and  the  flat  layer  theory 
to  yield  teleseismic  body  wave  radiation  from  explosion. 

Consider  plane  harmonic  P  waves  of  unit  amplitude  incident  at  an  angle  6  at  the  base  of  flat  layered 
media.  In  each  layer  m,  the  resulting  wave  field  will  be  a  superposition  of  up-  and  down-going  plane 
waves  and  corresponds  to  compressional  and  rotational  displacement  potentials  of  the  form 

(l)„,(x,z,to)  =  [A,mexp(-tkYamZ)  +  A2mexp(/kYamZ)]exp(tk(ct-t-x)] 

\|/m(x,z,(0)  =  [Bimexp(-ikY|3mZ)  +  B2raexp(ikYpmZ)]exp(tk(ct-i-x)]  (4) 

with 

k  =  (O/c 

Yam  =  (cW-l)‘'' 
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where  a  is  the  compressionai  velocity,  p  the  shear  velocity  and  c  the  phase  velocity. 

The  coefficients  Akm  and  Ekm  are  easily  determined  using  the  Thomson-Haskell  matrix  formulation 
(Thomson,  1950  ;  Haskell,  1953). 

If  u  denotes  the  resulting  horizontal  displacement  at  a  point  P  (Fig.  2)  the  reciprocity  theorem  implies 
that  a  horizontal  unit  force  applied  at  P  will  produce  at  Q  a  displacement  u  in  the  0-direction.  Differenti¬ 
ating  u  with  respect  to  x  will  then  yield  the  contribution  from  a  horizontal  unit  force  dipole  acting  at  P. 
Similarly,  the  contribution  from  a  vertical  force  dipole  will  be  given  by  the  z-derivative  of  the  vertical 
displacement  v  at  P.  Therefore,  a  spherical  explosive  source  with~spectral  seismic  moment  Mo((0)  will 
produce  at  Q  a  radial  spectral  displacement 


Ur(tO)  = 


3u 

L  ^ 


I]  Mo(«) 


(5) 


that  is,  X  and  z  being  the  source  coordinates 

2 

Ur(to)  =  Mo(a>)  [Ainiexp(-ikYo™z)  +  A2mexp(ikYamZ)]exp(/k(ct+x)]  (6) 


The  tangential  displacement  ue  can  be  obtained  in  the  same  way  by  computing  the  amount  of  dilatation 
produced  by  a  unit  force  applied  at  Q  along  the  direction  perpendicular  to  r .  The  incident  wave  is  now  an 
SV  wave. 

5.  Numerical  computations 

We  will  only  give  numerical  results  for  the  explosive  point  source  as  the  simplest  type.  But  they  will 
serve  as  guide  to  prove  how  to  simplify  the  forward  problem  using  the  reciprocity  theorem. 

For  a  teleseismic  body  wave  calculation,  an  explosive  point  source  is  set  in  a  layer  over  half-space. 
This  source  is  driven  with  a  Gaussian  time  function  given  in  Figure  3.  The  model  used  in  this  calculation 
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is  given  in  Table  1.  We  wish  to  compute  far-field  displacements  from  an  explosive  point  source  in  a  layer 
using  the  reciprocity  theorem.  In  the  reciprocal  approach,  plane  P  or  SV  wave  is  directed  from  the  half¬ 
space.  Using  the  Thomson-Haskell  matrix  method,  the  horizontal  and  vertical  transfer  functions  are 
computed  at  the  source  position  (Fig.  4  (a)  &  (b))  and  then  the  resultant  P  or  SV  wave  responses  are 
computed  by  applying  the  reciprocity  theorem  to  these  transfer  functions  (Fig.  4  (c)).  These  resultant 
responses  are  the  responses  at  the  half-space  caused  by  the  explosive  point  source  in  a  layer.  In  Figure  5, 
we  compare  P  wave  response  computed  from  the  reciprocal  approach  with  that  computed  from  the  for¬ 
ward  approach.  For  the  forward  problem,  we  used  Fuchs’  expressions.  He(1966)  derived  analytical 
expressions  for  the  transfer  function  for  dilatational  body  wave  radiating  into  the  mantle  from  a  point 
source  in  a  layered  crust  using  the  Thomson-Haskell  matrix.  This  result  shows  that  the  response  computed 
from  the  reciprocal  approach  is  equivalent  to  that  computed  from  the  forward  approach. 

On  the  basis  of  this  reciprocal  approach,  we  consider  the  influence  of  the  angle  of  incidence  into 
the  half-space  and  the  source  depth.  The  P  wave  responses  have  been  determined  for  angles  of  incidence 
ranging  from  5°  to  50°  at  a  fixed  source  depth,  2  Km  (Fig.  6).  It  may  be  concluded  that  the  P  wave  incident 
into  the  half-space  at  various  angles  does  not  vary  significantly  with  angle  of  incidence.  We  will  now 
place  the  source  at  six  different  depths :  0.5,  1,  2, 4,  7  and  10  Km  (Fig.7).  The  angle  of  incidence  into  the 
half-space  is  chosen  as  15°.  This  result  shows  that  the  P  wave  response  for  an  explosive  source  m  a  layer 
may  be  rather  sensitive  to  changes  in  the  depth  of  source. 

6.  Future  research 

(1)  The  reciprocal  approach  for  dislocation  source  in  homogeneous  layered  media 

(2)  The  reciprocal  approach  for  explosion  source  in  three-dimensional  heterogeneous  media 

(3)  The  reciprocal  approach  for  dislocation  source  in  three-dimensional  heterogeneous  media 

(4)  The  application  of  this  reciprocal  method  to  teleseismic  data 
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Forward  problem  Reciprocal  problem 


Figure  1 .  Reciprocity  principle.  Relation  between  the  force  f(t) 
and  an  induced  displacement  u(t). 
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Figure  2.  Geometry  of  the  problem 


20 


P-wave  velocity  S-wave  velocity  density  layer  thickness 
(Km/sec)  (Km/sec)  (g/cnP)  (Km) 


Table  1 .  Layer  over  half-space  model 


Time(sec) 


Figure  3.  Incident  wave  time  history. 
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DISPLACEMENTS  FOR  THE  EXPLOSION  SOURCE 
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Figure  4.  Horizontal  and  vertical  transfer  functions  computed  from  the  reciprocal  ap¬ 
proach  and  the  resultant  P  and  S  V  wave  displacements  computed  by  applying 
the  reciprocity  theorem  to  these  transfer  functions.  The  angle  of  incidence  is 
15°  and  the  source  depth  is  1  Km.  The  model  used  is  given  in  Table  1. 

(a)  the  transfer  functions  for  the  incident  P  wave. 

(b)  the  transfer  functions  for  the  incident  SV  wave. 

(c)  the  resultant  P  and  SV  wave  displacements. 
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P  WAVE  DISPLACEMENTS  FOR  THE  EXPLOSION  SOURCE 
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Figure  5.  A  comparison  of  the  synthetic  P  waves  computed  from  the  forward  problem 
and  from  the  reciprocal  problem.  The  angle  of  incidence  is  15  and  the  source 
depth  2  Km.  The  model  used  is  given  in  Table  1. 
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Time(sec) 


Figure  6.  P  wave  responses  at  the  source  depth,  2  Km,  at  various  incidence  angles.  The 
angles  of  incidence  range  from  5°  to  50°  with  interval  5°.  The  model  used  is 

given  in  Table  1. 


Source  depth  =  0  Km 


Figure  7.  P  wave  responses  at  the  angle  of  incidence,  15°,  at  various  source  depths.  The 
source  depths  are  0 , 0.5,  1,  2»  4,  7  and  10  Km.  The  model  used  is  given  in 
Table  1. 
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We  investigate  theoretically  the  effects  of  anisotropy  by  calculating  synthetic  seismograms  for 
simple  plane-layered  structures  with  a  free  surface.  For  the  theoretical  part  of  our  work,  we  take  as  our 
basis  the  work  of  Taylor  &  Crampin  (1978). 

(1)  Theoretical  aspects 

We  consider  an  n  layered  half-space,  with  the  axes,  and  the  layer  and  interface  numbering  arranged 
as  Fig.  1.  The  direction  of  the  apparent  velocity  c  along  the  surface  is  confined  to  the  Xi  direction,  and 
all  the  elastic  tensors  are  rotated  accordingly. 

From  the  elastodynamic  equations  of  motion,we  obtain  a  matrix  equation  for  p,  frequently  called  the 
slowness  equation: 

Fa  =  (Rp^  4-  Sp  +  T  -pc“I)a  =  0,  (1) 

Since  R  is  non-singular,  equation  (1)  can  be  written  in  partitioned  matrices  as  a  linear  eigenvalue 
problem  for  p: 


-R-‘S 

I 


(2) 


The  solution  of  the  slowness  equation  in  the  form  of  the  eigenvalue  problem  and  simple  scaling  of 
the  eigenvectors  using  the  orthogonality  of  the  row  and  column  vectors  for  eigenvalues  allows  us  to 
construct  the  propagator  matrices  of  Gilbert  &  Backus  (1966)  without  numerical  inversion. 

The  product  of  the  propagator  matrices  relates  the  displacement-stress  vectors  at  the  top  of  the 
half-space  to  those  at  any  other  interface.  The  solution  is  determined  by  relating  the  half-space  vector 
to  the  boundary  conditions  at  the  free  surface,  where  the  components  of  normal  stress  vanish.  We  have 

(6,k.  52k,  63k,/4,/5,/6)n^  =  En''  G  (U,,  Uj,  U3,  0,  0,  0)o''  (3) 

where  G  =  O'  m=n-l  Gfm 

(2)  Results  of  calculations 

The  incident  pulse  to  be  transmitted  through  the  structure  from  the  lower  half-space  is  given  by 
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Fig.  2.  In  this  paper  the  models  given  by  Table  1  are  investigated.  The  anisotropic  media  in  these 
models  are  transversely  isotropic  medium  (Table  2  (a))  and  the  anisotropic  medium  with  the  (OOl)-cut 
orthorhombic  olivine  (Table  2  (b)).  In  our  models  we  take  the  axis  of  symmetry  in  the  plane  (xi,  X2), 
making  an  arbitrary  angle  (j)  with  the  xi-axis,  i.e.,  to  the  horizontal  projection  of  the  wave  propagation 
direction.  The  angle  (|)  is  reckoned  from  the  xpaxis  in  an  counter-clockwise  sense  when  viewed  from 
above. 

Figure  3  shows  the  dependence  of  the  velocity  of  each  body  wave  on  the  propagation  direction  in 
the  transversely  isotropic  medium  and  the  anisotropic  medium  with  (OOl)-cut  orthorhombic  olivine 
given  in  Table  2. 

We  now  examine  the  synthetic  seismograms  computed  for-the  models  in  Table  2  with  plane  body 
waves  of  the  three  types  P,  SV,  and  SH  approaching  from  the  isotropic  half-space.  The  most  conunon 
property  of  all  in  the  synthetic  seismograms  is  the  presence  of  the  traces  of  the  three  wave  types,  P,  Si, 
and  $2  on  every  component  (Fig.  4  &  5).  While  propagating  along  the  layer,  the  P,  Si,  S2  waves,  by 
virtue  of  their  velocity  differences,  arrive  at  different  times  and  manifest  themselves  as  separate  inci¬ 
dences  on  the  seismograms.  For  the  case  of  the  incident  P-wave,  the  presence  of  anisotropy  couples 
the  P,  SV  and  SH  wave  motion  so  that  P  waves  incident  on  the  anisotropy  layer  from  below  produce 
P,  SV  and  small-amplitude  SH  waves  at  the  surface.  In  the  case  of  incidence  of  the  SV  and  SH  waves 
on  the  lower  interface,  the  anisotropy  of  the  layer  manifests  itself  most  distinctly  in  the  splitting  up  of 
the  SV  and  SH  waves  into  S|  and  S2  waves  at  the  lower  interface,  the  latter  having  different  velocities. 
The  ratio  of  the  amplitudes  of  the  S|  and  S2  waves  strongly  depends  on  the  angle  ((). 

From  these  results,  we  can  infer  that  the  effects  of  anisotropy  on  incident  P  waves  are  small,  while 
the  effects  on  incident  S  waves  are  comparatively  large. 
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Notations 


a  the  amplitude  vector,  with  elements  { aj } ,  of  a  particular  plane  wave  decomposition 

Cjkmn  the  elements  of  the  elastic  tensor,  not  necessarily  referred  to  the  principal  axes 

/  the  vector  of  excitation  functions  with  elements  {^},  j=l,2,....,6 


I 

P 

R,  V,  and  T 

S 

T" 

5jm 

P 

CO 

G 


the  3  X  3  identity  matrix 
the  normalized  slowness  vector 

submatrices  of  the  full  elastic  tensor,  with  elements  {cjska},  {cjiks},  and  {cjiki }, 
respectively 

=  V  +  V'^ 

=  T-pch 


the  Kronecker  delta  function 
density 

the  angular  frequency 

the  propagator  matrix  which  is  expressed  by  Gn  =  EnDnEn ' 
where 


E 


n 


(0  i\(ap  a'p'A 
\R  y)\  A  A'  J 


T 


Dn  =  diag[exp(-/to/7jdn/c)] 


dn  the  thickness  of  the  nth  layer 

P  =diag(/7i,;?2,P3) 


P'  =  diag(/?4,  P5,  Pe) 

A  =  (ai,a2,a3) 


A'  =  (a4,a5,a6) 
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Table  1.  Structural  models  used  for  the  calculations. 


Model  No. 

Layer  No. 

Vp 

(km/sec) 

Vs 

(km/sec) 

Density 

(gm/cc) 

Thickness 

(km) 

1 

1 

transversely  isotropic  medium 
(Table  2) 

3.3 

35.0 

2 

8.12 

4.7 

3.3 

half-space 

1 

7.25 

4.2 

2.9 

30.0 

2 

2 

((X)l)-cut  olivine  (Table  2) 

3.324 

30.0 

3 

10.0 

5.7 

3.6 

half'Space 

Table  2.  (a)  Elastic  constants  of  a  transversely  isotropic  medium. 

A  =196.0  C  =  244.0  F  =  70.0  L  =  82.5  N  =  67.0  x  lO’ N-m'^ 
(b)  Elastic  constants  of  onhorhombic  olivine. 


ciiii  =  324.0 

C2222  =  198.0 

C3333  =  249.0  X  10^  N-m 

Cl  122  =  59.0 

C2233  =  78.0 

C33n  =  79.0 

Ci2i2  =  79.3 

C2323  =  66.7 

C1313  =81.0 
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Figure  1.  Coordinate  axes  for  multilayered  models  where  the  layers  can  be  isotropic  or  anisotropic. 
0  is  the  angle  of  incidence. 
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Figure  4.  Displacements  of  the  free  surface  produced  by  (a)  P,  (b)  SV,  (c)  SH  at  25®  angle  of  incidence 
on  transversely  isotropic  layer  in  Model  1  of  Table  1.  (j)  is  angle  between  the  axis  of  symmetry 
of  the  transversely  isotropic  layer  and  the  horizontal  projection  of  the  wave  propagation  direction. 
The  input  pulse  used  here  is  given  by  figure  1(a). 
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Displacements  of  the  free  surface  produced  by  (a)  P,  (b)  SV,  (c)  SH  at  35°  angle  of  incidence  on  (OOl)-cut  olivine  in  Model  2 
of  Table  1.  (|)  is  the  angle  measured  from  (100)  towards  (010).  The  input  pulse  used  here  is  given  by  figure  1(b). 
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Introduction 

The  purpose  of  this  tutorial  is  to  provide  some  documentation  of  the 
codes  in  the  TELEDB  system  so  a  reasonably  proficient  seismology 
graduate  smdent  can  quickly  perform  moment  tensor  inversions  using  IRIS 
waveform  data.  I  had  high  hopes  of  integrating  my  synthetic  waveform  md 
inversion  codes  with  the  JSPC  software  package  DATASCOPE,  a  fine 
database  system  written  by  Dan  Quinlan  at  the  JSPC.  As  it  turned  out, 
because  of  the  amount  of  work  that  I  put  into  rewriting  the  teleseismic 
waveform  program,  TELEDB  falls  far  short  of  utilizing  the  full  range  of 
utilities  and  formats  that  is  in  DATASCOPE.  The  current  version  of  Ae 
software  incorporates  some  very  useful  DATASCOPE  utilities  for  building 
database  files  from  IRIS  SAC  files.  I  tried  to  preserve  the  integrity  of  the 
CSS3.0  database  schema,  but  as  time  pressed  for  completing  programming 
for  this  short  course,  I  had  to  "bend"  it  quite  a  bit  by  making  up  new  file 
formats  to  interface  with  various  waveform  utilities,  the  synthetics  code,  and 
the  inversion  code.  Although  I  would  like  to  more  fully  integrate  this 
package  with  DATASCOPE,  perhaps  some  enterprising  graduate  student 
might  be  better  inclined  and  more  proficient  at  programming  than  me  to  do 

so.  ,  ^  • 

Aside  from  several  DATASCOPE  modules,  my  style  of  programming 

consists  of  patching  together  several  FORTRAN  77  programs  through  a 
UNIX  C  shell  script.  Most  scripts  set  up  a  small  parameter  file  that  is  then 
used  by  one  or  more  fortran  programs  that  then  access  database  files  and 
waveform  files.  Use  of  SAC  (Seismic  Analysis  Code,  Lawrence  Livermore 
Laboratory)  occurs  often  in  these  scripts  and  the  data  are  assumed  to  reside 
in  SAC  binary  files.  A  particular  directory  structure  is  assumed  by  these 
scripts  and  must  be  set  up  by  the  user  using  a  utility  script.  I  am  sure  there 
are  more  elegant  ways  of  programming,  but  I  have  always  l^en  attracted  to 
the  simple  and  the  direct,  knowing  how  easy  it  is  to  make  mistakes. 

Before  getting  into  the  details  of  the  software  system,  the  next  section 
outlines  theory  for  computations  of  the  Green's  functions  for  teleseismic 
body  waves  and  the  waveform  inversion. 

Computation  of  Teleseismic  Body  Waves 

Accurate  moment  tensor  inversion  of  seismic  waveforms  requires  that 
the  Green's  functions  be  known.  This  is  actually  the  hardest  part  of  all 
seismological  problems  since  knowledge  of  earth  strucmre  is  often  very 
sketchy.  Fortunately,  teleseismic  body  waves  can  be  explained  with  fairly 
simple  theory.  Rays  from  a  buried  source  interact  with  structure  near  the 
source  but  then  are  largely  unaffected  by  propagation  through  the  lower 
mantle  until  arriving  at  the  receiver.  Receiver  effects  for  vertical  P  waves 
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and  tangential  S  waves  are  relatively  minor  and  do  little  to  mask  the  signal 
controlled  mostly  by  source  radiation  pattern  and  near-source  structure. 

A  FORTRAN  77  program  called  dishask  ( I  will  use  bold  characters 
for  modules  in  TELEDB)  performs  the  computations  for  the  problem 
outlined  in  Figure  1.  The  far-field  response  of  a  point  source  m  plane 
layered  media  is  calculated  using  the  Thomson-Haskell  technique.  The 
phase  term  for  propagation  and  geometrical  spreading  in  the  halfspace  is 
replaced  by  the  geometrical  spreading  effect  of  a  P  or  S  wave  propagating 
through  a  reference  spherical  earth  structure.  The  far-field  source  response 
is  then  input  at  the  receiver  side  as  a  plane  wave  incident  under  a  plane 
layered  receiver  structure.  The  program  allows  for  differing  source  and 
receiver  structures.  The  source  structure  may  include  a  surface  water  layer 
for  modeling  events  under  oceans.  Receiver  structures  are  allowed  to  be 
different  since  stations  are  situated  in  many  different  tectonic  environments 
and  typical  structures  can  have  significant  effects  on_the  resulting  response, 
depending  on  displacement  component.  The  TELEDB  database  structure  is 
setup  so  that  as  knowledge  of  structure  under  various  receivers  becomes 
available,  these  structure  models  can  then  be  incorporated  in  future 
o^IcuIh  lions 

FoUowing  the  notation  of  Harkrider  (1964;  1970)  the  downgoing  P 
and  SV  wave  potential  coefficients  are  given  by  the  following  matrix 
equations.  Geometry  for  the  source  problem  is  shown  in  Figure  2. 
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Figure  1 
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A„  =a„  a„  •••a„  (4) 

^Sl  ^Sl  ^S-l  ‘S 

60  -  source  jump  conditions 


and  aRi  is  the  Haskell  layer  matrix. 
At  the  free  surface 


cr„  =0  . 


(5) 


For  a  water  layer  on  top  of  the  layered  stack  with  a  thickness  of 


(6) 

c  -'ll  c 


a  =A  (7) 

^0  J  21  C 


Af  is  the  Haskell  layer  matrix  for  a  fluid. 

The  equations  for  the  downgoing  SH  wave  potential  are  given  by 


/ 

n 

f 

nj 


=J, 


0 


+  A‘ 


i^L  ^ 

6 

I  "  y 

Str 


(8) 


where  Jl  and  Al'^  are  defined  in  an  analogous  way  to  the  P-SV  system 
above. 


The  source  jump  conditions  are  given  by: 


^  ; 


43 


TELEDB  Tutorial  Page  7 


The  source  coefficients  etc.,  are  given  by 
Vertical  Strike-Slip:  v 

a 

4  =-'5^=-^* 
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^ 

O  —  ^ 

•^03  “  "^03  “  r 

(20) 

Vertical  Dip-Slip 

*^01  “  “*^01  “ 

(21) 

^  KKkl  -  Ik^ ) 

^02  -  ^02 

(22) 

So3  =  “-^03  “  ~^^fi 

(23) 

Compensated  Linear  Vector  Dipole  (CLVD) 

^  Ki(2k:l-3k^ 

*^01  “  *^01  f 

a 

-  (24) 

So2  =  “■^02  “ 

(25) 

■^03  "  “^03  “  ^ 

(26) 

Isotropic  Source  (Implosion) 

,  -  _  Kico^ 

*^01  “  ‘^01  ~  2 

«  ''a 

(27) 

*^02  “■^02“^ 

(28) 

(29) 

where. 
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K  = 


flLHD(co) 

AnpaP' 


(30) 


ik„=-  (31) 

«  a 


k  =- 

p 


(32) 


Spectral  displacements  in  the  halfspace  are  given  by 


-ik„  R 


a 

n 


(33) 


Ir.  ,,  -V 

^SV  k  ^  R 


(34) 


-ik.  R 

Pn 


^SH  =  — 


(35) 


1=2  for  vertical  strike-slip 
=  1  for  vertical  dip-slip 
=  0  for  CLVD  and  isotropic  sources. 
Solution  of  the  matrix  equations  give  the  downgoing  wave  coefficients: 

=^-L-[A  +  B  +  C  +  D]  (36) 

"  detj^ 

m'  =—L-[E  +  F  +  G  +  K]  (37) 

"  detj^ 
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H,  J,  -H.J, 

A'  ni  ^1 
"  r  -  y 


(38) 


where 


-1 


“  ^R^R, 


SI 


(  Ur  ^ 

6 

I 


V  y 


Hi,=JLA2^, 


5(7, 


Sr, 


■  (^L  ^ 

5 

[  ^  J 


Sr, 


(39) 


(40) 


B  =  {(■/(!„  -  )-^R,3  -  (•^R„  -  \ 


C  =  H 


(43) 

(44) 

(45) 


47 


TELEDB  Tutorial 


Page  11 


All  of  these  quantities  can  be  explicitely  found  in  the.  dishask  source  code. 

Geometrical  spreading  is  calculated  for  a  flattened  spherical  earth 
model  using  the  following  equation. 


2  -  1  P  ( 


^0  = 


\  dbcAsin0. 


(50) 


where. 


Pj  =  density  of  mantle  at  source 
Pq  =  density  of  mantle  at  receiver 
Vq  =  wave  velocity  at  source 
Vj  =  mantle  wave  velocity  at  receiver 

X  -  distance 
p  =  ray  parameter 

9  =  spherical  earth  distance  (radians) 

%  • 

‘  n- 

ao  is  the  amplitude  of  the  P  or  S  wave  in  the  receiver  inantle  directly  below 
the  layered  structure.  There  is  a  series  of  subroutines  in  dishask  to 
determine  the  appropriate  ray  parameter  and  geometrical  spreading  given  a 
source-receiver  distance. 
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Finally,  the  receiver  effect  is  computed  using  the  Thomson-Haskell 
method  for  an  incident  plane  wave  under  a  plane  layered  structure.  The  free 
surface  displacements  are  given  by: 


detj^ 

_  2c 

detj^ 

c 

detj^ 

c 

detj^ 

SH  A  .  A 

p  B  r„  A.  +Aj 


(51) 

(52) 

(53) 

(54) 

(55) 


where,  w,u,v  sltq  vertical,  radial  and  tangential  displacements, 
respectively.  The  notation  for  the  layered  receiver  structure  follows  the 
same  conventions  as  the  source  problem  and  Figure  2,  except,  obviously, 
there  are  no  source  layers. 

Complex  spectra  for  the  4  Green's  functions  are  computed  and 
multiplied  with  the  spectra  for  the  appropriate  receiver  response  and  the 
geometrical  spreading  factor.  An  inverse  FFT  brings  the  spectra  into  the 
time  domain. 

There  is  one  additional  detail  which  may  be  useful  to  researchers. 
Attenuation  is  included  in  both  source  and  receiver  structures  through  the 
use  of  complex  wave  velocity.  "Effective  velocity",  veff,  is  given  by: 


49 


TELEDB  Tutorial  Page  13 


where  v  is  P  or  S  wave  velocity  and  Q,  P  or  S  wave  Q.  coq  is  an  arbitrary 
attenuation  band  cutoff  frequency  (set  here  to  a  period  of  1000  sec)  and  y, 
Euler’s  number.  This  form  of  complex  velocity  uses  Futterman's  (1962) 
result  for  phase.  Attenuation  due  to  wave  propagation  in  the  mantle  also 
uses  the  Futterman  operator  and  is  included  in  the  final  spectral 
multiplication. 

The  source  time  function  is  parameterized  by  a  series  of  triangular 
elements  and  is  discussed  in  the  next  section.  I  originally  used  boxcar 
elements  in  past  inversions  but  Nabelek’s  (1982)  triangular  element 
parameterization  is  more  appealing  since  it  produces  a  smoother  time 
function. 


Inversion  of  Teleseismic  Body  Waves 


The  inversion  for  moment  tensor  and  time  function  follows 
Wiggin's(1972)  method  of  Singular  Value  Decomposition  and  singular  value 
truncation  for  damping.  Because  joint  inversion  of  both  sets  of  parameters  is 
nonlinear,  I  discuss  the  details  of  the  source  parameterization  in  addition  to 
those  in  the  shortcourse  tutorial  on  moment  tensors. 

Displacements  can  be  written  in  terms  of  the  moment  tensor 


(57) 


where 


(58) 


are  the  moment  tensor  elements,  and  G.. .  are  the  Green’s  functions.  In  a 
cylindrical  coordinate  system,  and  for  a  deviatoric  moment  tensor  source 


(t,r,z)A'.i6,X,8)  (59) 
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1  <  * 


w{t,r,z,Az)-s{t)*^H'^  (t,r,z)M'.  (62) 

.  -I  * 


where 


and 


Ml  Mu 

"2=^22 

M'=M,2 

a/;=m,3 

"5-"23 


(63) 


H'  =--//  cos(2Az)+^/^ 

W,  2  2  ^2 

//'  =-//  cos(2 Az)+—H 

2  2  ^3 

H'  =-H  sin(2Az)  (64) 

W3  w, 

H'  =H  cos(Az) 

H'  =H  sin(Az) 

Wj 

The  Green's  functions,  e.g.,  Hy^iy  are  given  above  for  a  vertical  strike-slip 
source  (i=l),  vertical  dip-slip  source  (i=2),  and  CLVD  source  (i=3). 

The  source  time  function  parameterized  by 

(65) 

ifc=l 
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where  B(t  -  Tk)  is  a  triangular  element  characterized  by  equal  rise  times  and 
fall-off  times  and  a  width  of  Ax. 


Tk  =  —AT{k  - 1 )  (66) 

The  far-field  time  function  is  normalized  such  that 

Sk  =  I  (67) 

2  k=I 


Then, 

L  5  ^  ~ 

w(r,r,z,Az)=[2^5^B(/- 

ife=l  i=l  ' 

or 


w(/,r,z,Az)=x|;5^M;[B(r-  (r,r,z,Az)]  (69) 

ik=ij=i 

The  object  of  the  inversion  is  to  find  the  moment  tensor  and  time 
function  coefficients.  Since  the  equation  shows  a  direct  tradeoff  between 
moment  tensor  and  time  function,  an  equation  of  constraint  is  added  to  the 
final  matrix  equations  before  inversion.  This  equation  states  that  the  time 
function  perturbations  average  to  zero  so  that  the  time  function  remains 
normalized: 


Asife  =  0  (70) 

k=l 

The  matrix  equation  are  set  up  in  a  standard  way 


where 
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O .  =  observed  ground  motion  at  time  j 

C®  =  starting  model  synthetic  at  time  j  (72) 
AP.  =  parameter  perturbation 


IpO 


Weighting  matrices  are  also  included  where,  W  is  the  parameter  weighting 
matrix  and  S  is  the  data  covariance  matrix 

S..=cov(0;0.]  (73) 


I  assume 


S„  =  af  (74) 

The  transformed  system  of  equations  become 

a'ap'  =  Ac'  (75) 

where 

a' 

Ap'=W"^/^AP  (76) 

AC'=S"^/^AC  . 
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Using  Singular  Value  Decomposition 

a'  =  UAV^  (77) 

AP'  =  VA^'U^'ac’  (78) 

The  variance  of  the  solution  is  examined  to  determine  if  Aere  are  any 
near-zero  singular  values  which  would  indicate  and  unstable  inverse.  An 
arbitrary  cut-off  threshold  is  specified  to  cut  off  these  small  singular  values 
which  damps  the  inverse  solution.  Since  small  singular  values  are  usually 
associated  with  the  time  function  parameterization,  another  equation  of 
constraint  can  be  included  in  the  system  of  equations  prior  to  inversion, 
which  specifies  that  everyother  time  function  element  is  inverted  with  the 
element  between  being  specified  by  the  average  of  the  adjacent  elements. 
This  averaging  technique  reduces  the  number  of  independent  variables  and 
helps  convergence  sometimes.  Note  that  you  must  specify  an  odd  number  of 
time  function  elements  to  perform  the  averaging. 

Singular  value  decomposition  allows  a  number  of  useful  secondary 
results  to  be  calculated  during  the  inversion.  Data  importence  is  evaluated 
by  summing  along  the  columns  of  the  information  density  matrix  for  the 
seismogram  of  a  particular  station.  The  information  density  matrix,  or 
resolution  matrix  is  given  by  U.  The  parameter  resolution  matrix  V  V 
is  also  useful  to  examine  for  the  tradeoff  between  parameters.  If  all  singular 
values  are  kept,  this  matrix  is  the  identity  matrix  and  all  parameters  are 
perfectly  (in  linear  theory)  resolved,  aside  from  errors  induced  by  noise  in 
the  data. 

The  success  of  the  inversion  is  examined  several  ways.  One  is  to 
simply  see  how  good  the  data  waveforms  have  been  fit.  This  measure  is 
calculated  by; 


(79) 

The  least-squares  error  is  also  calculated  to  examine  how  good  the 
inverse  approximated  the  data  perturbations: 

lAAP  -  AC  I  (80) 

In  addition  to  these  formal  error  criteria,  the  source  model  is 
decomposed  into  the  major  double  couple  and  a  remainder  CLVD  source. 
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The  amount  of  non-double  couple  source  can  be  interpreted  as  a  measure  of 
misfit,  if  you  believe  that  earthquake  sources  must  conform  to  the  double 
couple  or  dislocation  model. 

In  practice,  it  is  efficient  to  first  perform  a  linear  moment  tensor 
inversion,  assuming  a  source  time  function  and  source  depth,  to  get  a 
reasonable  starting  model  for  the  source  radiation  pattern.  Once  an 
approximate  moment  tensor  has  been  found,  joint  inversion  for  time  function 
and  moment  tensor  is  performed  since  this  is  a  non-linear  process.  Sever^ 
interations  may  be  needed  for  convergence,  if  it  occurs.  The  l^st  depth  is 
chosen  by  examining  the  three  error  criteria  and  minimizing  all,  if  possible. 

The  TELEDB  System 


Overview 

The  TELEDB  software  tools  have  been  designed  for  the  formidable 
problem  of  downloading  large  amounts  of  IRIS  waveform  data  from  the 
SPYDER  system  in  SAC  binary  format  ,  preparing  the  data  for  moment 
tensor  inversion,  and  then  inverting  the  data.  The  waveform  data  are 
interogated  using  DATASCOPE  software  to  produce  database  files,  and 
these  files  are  used  in  subsequent  processing  steps. 

The  large  number  of  waveform  files  available  for  any  one  event  is  a 
boon  for  science  but  makes  the  data  processing  extremely  tedious  unless  big 
pieces  of  the  processing  can  be  automated.  The  approach  used  here  is  to 
make  use  of  database  files  to  keep  track  of  the  data  and  data  processing  and 
to  standardize  the  file  system  so  that  scripts  can  be  used  over  and  over  for 
different  earthquake  events. 


File  Structure 

A  basic  UNIX  directory  structure  is  assumed.  Assume  that  you  are 
working  on  an  earthquake  from  the  Yunnan  province  in  China.  The  root 
event  directory  may  be  called  "yunnan"  with  various  subdirectories  created 
underneath  this  root. 

./yunnan  root  directory 

./yunnan/data  root  data  directory 

./yunnan/data/raw  raw  data 
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./yunnan/data/process 

processed 
data  of  several 
types 

./yunnan/data/invert 

final 

windowed  and 
time  shifted 
waveform  data 

./yunnan/sacplz 

SAC  polezero  files 

./yunnan/synth 

An  extra  directory 
for  placing 
synthetic 
seismograms 

./yunnan/greens 

inversion  Green's 
functions 

./yunnan/invert 

inversion  results 

./yunnan/db 

database  files 

./yunnan  may  be  a  subdirectory  under  "teledb"  which  contains  the 
executables,  source  code,  earth  models,  and  SEED  response  files. 

teledb 


teledb/bin  Shell  scripts  and 

binaries 


teledb/src 

teledb/emodels 

teledb/responses 


source  code 

earth  model  files 

SEED  instmment 
response  files 
downloaded  from 
IRIS 


Most  all  scripts  require  that  a  UNIX  environmental  variable  called 
$PSUMOM  be  set  in  your  login  or  .cshrc  file.  For  example,  if  the  teledb 
directory  is  under  /tremor/sO/cal,  then: 
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setenv  PSUMOM  /tremor/sO/cal/teledb 

After  setting  PSUMOM,  an  event  directory  with  subdirectories  can  be  built 
by  invoking  the  C  Shell  Create_directory.csh  which  is  self-explanatory. 
(See  the  Appendix  for  C  Shell  script  listings.) 

Making  a  Database 

After  the  basic  directory  structure  has  been  set  up,  obtain  data  from 
IRIS.  The  system  was  developed  around  the  practical  problem  of 
downloading  data  from  the  SPYDER  system  over  INTERNET.  Log  into  the 
SPYDER  system,  choose  an  event,  and  download  SAC  binary  data  into  the 
/data/raw  subdirectory. 

After  the  data  have  been  obtained,  make  sure  there  are  no  other  files  in 
the  "/data/raw"  directory  so  database  files  may  be  built  using  the  sac2db 
command.  If  you  want  database  tables  named  with  the  root  "yunnan",  for 
example,  you  would  type,  in  the  directory: 

sac2db  *  yunnan 

This  program  finds  all  the  files  in  the  directory  and  builds  a  CSS3.0 
database.  Move  the  database  files  to  the  7db’’  directory.  At  this  point  you 
will  need  to  make  a  good  origin  file,  since  sac2db  will  have  made  one  that 
contains  a  lot  of  garbage  because  IRIS  SAC  files  do  not  contain  all  of  the 
necessary  information.  This  can  be  done  using  the  DATASCOPE  utility 
dbe.  In  the  "/db"  directory  type,  for  example: 

dbe  yunnan 

Get  up  the  menu  and  make  a  new  origin  table.  Save  it. 

Instrument  Corrections 

Now  that  the  necessary  database  tables  have  been  made,  notably, 
"event.wfdisc"  and  "event.origin",  the  data  can  be  processed  for  inversion. 
This  consists  of  the  following  steps: 

1.  Calculate  P  and  S  wave  arrival  times  for  phase 

interpretation.  Put  these  times  into  the  data  file. 

2.  Constmct  a  SAC  polezero  file  for  a  particular  data 

trace. 

3.  Write  a  SAC  macro  file  for  instrument  correction. 
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4.  Correct  the  data  for  instrument  response  and  make 
ground  displacements. 

These  functions  are  contained  in  the  C  Shell  prewave.csh.  A  reading 
of  this  script  shows  that  it  builds  an  input  file  for  a  Fortran  77  program 
called  prewave  that  interogates  the  database  files  and  SAC  data  files  for 
appropriate  information,  prewave  creates  another  C  Shell  which  generates 
SAC  polezero  files  from  the  SEED  responses  contained  in  teledb/responses 
using  the  program  evalresp,  which  was  obtained  from  the  IRIS  f^  archives. 
Note  that  prewave  looks  for  an  appropriate  response  file,  but^  if  it  caimot 
find  one,  it  will  discard  the  particular  data  waveform  from  consideration.  It 
is  the  job  of  the  analyst  to  make  sure  that  SEED  response  files  for  each 
waveform  channel  are  available  in  the  "teledb/responses"  directory.  These 
are  available  from  a  SEED  tape  volume  if  the  data  are  obtained  that  way.  I 
am  providing  the  most  recent  set  of  responses  available  from  the  IRIS  ftp 
Qrcliivcs. 

AU  of  the  C  Shell  scripts  are  fairly  self-explanatory  if  they  are  invoked 

without  arguments.  ,1.^1 

The  prewave.csh  script  will  also  make  a  set  of  new  database  files, 

notably  newdb.wfdisc  which  contain  all  of  the  information  needed  on  data 
waveforms  that  have  been  corrected  for  instrument  response.  The  corrected 
data  will  also  be  written  to  the  "/data/process"  directory  for  further 
processing. 

Vector  Rotation 

The  next  step  in  processing  is  to  rotate  the  horizontal  waveforms  into 
the  theoretical  backazimuth  to  the  event,  to  produce  SV  and  SH  motions. 
This  is  done  using  rotate.csh.  This  C  Shell  builds  a  simple  input  file  for  the 
rotate  program,  rotate  then  interogates  the  database  .wfdisc  file  to  build  a 
SAC  macro  file  for  the  vector  rotation.  SAC  is  then  invoked  to  do  the  actual 
work.  New  database  files  are  produced  after  this  step,  also. 

Analyst  Decisions 

Until  now,  the  data  processing  has  been  routine  with  little  analyst 
interaction,  but  this  is  the  appropriate  point  to  really  examine  the  data  to  see 
if  there  is  anything  worthwhile  studying!  Use  SAC  and  look  and  P  and  SH 
wave  data.  How  is  the  signal-to-noise  ratio?  What  kind  of  high  pass  filter 
will  remove  baseline  and  long  period  noise  problems?  Which  stations  will 
give  useful  data?  Can  you  see  P  arrivals,  S  arrivals  based  on  the  travel  time 
estimates?  (P  and  S  times  will  be  marked  by  SAC  T1  and  T2  markers.) 
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Choose  your  data  wisely  and  jot  down  the  channels  and  stations  that  you 

want  to  forther  consider.  ,,1.1 

Another  C  Shell  script  is  invoked  to  create  the  next  database  table, 

which  I  call  a  ".stachan  table",  that  will  be  used  in  synthetic  seismogram 
calculations  and  in  the  moment  tensor  inversion.  This  is  where  I  diverge 
from  DATASCOPE  CSS3.0  conventions.  The  ".stachan"  file  can  be 
considered  to  be  a  parameter  file  in  the  DATASCOPE  scheme  of  things. 
The  C  Shell  is  called  create_stachan.csh  and  is  a  bit  tedious  since  it 
interactively  asks  you  questions  about  your  data.  Again,  it  constructs  a 
simple  data  file  which  is  read  by  a  fortran  program  called  sctable,  which 
actually  does  the  asking.  A  future  improvement  would  be  to  use 
DATASCOPE  itself  to  create  the  table,  but  that  entails  modifying  the 
CSS3.0  schema.  The  ".stachan"  file  is  a  basic  data  file  containing 
information  about  each  waveform  channel  of  interest. 


Data  Windowing  and  Decimation 

Next,  I  recommend  that  the  final  data  be  windowed  for  use  in  the 
inversion  program.  The  inversion  program  momvert  only  allows  data  array 
lengths  of  512  points,  so  the  waveform  data  and  synthetics  must  be  cut 
down,  and  often,  decimated.  Care  must  taken  to  make  sure  that  the  final 
data  seismogram  have  the  same  sampling  interval  as  its  corresponding 
synthetic  seismogram.  So,  at  this  point  you  must  carefully  determine  what 
sampling  interval  you  will  use  in  constructing  synthetics  and  how  to 
decimate  synthetic  and  observed  data  so  individual  data  channel/synthetic 
pairs  agree  in  sampling  interval.  This  can  be  a  problem  because  it  is  easy  to 
use  a  mixed  data  set  consisting  of  BH  (20  sps)  and  LH  (1  sps)  waveforms. 
The  utilities  for  dealing  with  this  problem  assume  that  you  may  BH  and  LH 
data  together.  You  are  on  your  own  with  other  data  sets. 

Suppose  you  will  be  constructing  synthetic  seismograms  with  a 
sampling  interval  of  0.2  sec,  or  5  sps.  That  means  you  will  have  to  deciniate 
the  observed  BH  waveform  data  using  a  decimation  factor  in  SAC  of  4. 
Conversely,  you  will  have  to  decimate  the  synthetic  seismograms  by  a  factor 
of  5  to  get  them  down  to  1  sps,  appropriate  for  the  LH  data  on  hand,  so  boA 
data  and  synthetics  may  need  decimation.  It  is  much  more  convenient,  the 
way  things  are  set  up  here,  to  calculate  the  synthetic  seismograms  with  a 
single  sampling  interval,  although  this  could  be  changed  with  a  little  work. 

A  little  more  analysis  with  SAC  is  necessary  before  windowing  and 
decimation  can  occur.  Bring  up  all  of  your  good  waveforms  in  a  SAC 
window  and  pick  the  onset  time  of  the  P  or  S  wave.  Be  careful  because  this 
is  a  critical  analysis  choice  that  will  affect  the  stability  and  veracity  of  the 
inversion.  Onset  times  are  marked  using  the  TO  marker  in  the  "ppk" 
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command  of  SAC.  Remember  to  write  over  your  data  file  to  save  the  time 

The  C  Shell  window_data.csh  is  invoked  to  perform  the  windowing. 
Your  choice  of  time  window  is  input  as  a  parameter  and  the  time  window  is 
constructed  to  start  5  seconds  before  the  TO  marker.  Decimation  factors  are 
also  input  as  well  as  a  "factor"  parameter  which  can  be  used  to  change  Ae 
units  of  the  data.  Most  GSN  instrument  responses  correct  to  m/sec  units, 
and  after  intergration,  to  m.  dishask,  the  synthetics  program,  uses  micron^ 
Thus,  a  factor  of  10^  would  be  usefully  multiplied  with  the  data  to  match 
units’.  (As  an  aside,  I  noticed  in  a  trial  of  this  procedure  that  not  aU 
instmment  responses  were  appropriately  normalized.  Some  data  amplitudes 
were  off  by  very  large  numbers.  Response  processing  should  be  carefully 

checked.} 

window_data.csh  wiU  create  a  SAC  macro  file  which  will  be  used  by 
SAC  to  process  the  data  yet  again.  These  final  data  waveforms  will  be 
•  stored  in  the  data/invert  subdirectory  for  use  with  momvert.  The  data 
waveform  files  will  have  the  suffix  ".win'  appended  to  the  original  name. 


Computation  of  the  Greeen's  Functions 

Create  synthetic  Green’s  functions  using  dishask.csh.  In  the  process 
of  using  the  create_stachan.csh  script,  you  will  have  t^en  asked  to  provide 
receiver  earth  structure  files  for  each  station.  A  few  sunple  earth  structure 
files  are  provided  under  "teledb/emodels/receiver  .  These  consist  of 
"shield",  "tectonic"  and  "oceanic"  models.  Each  model  file  has  the  suffix 
".res"  appended  to  it  denoting  receiver  g,arth  structure.  There  are  similar 
source  earth  models  in  teledb/emodels/source  except  the  files  have  the  suffix 
".ses"  for  source  g.arth  structure.  The  PREM  model  is  contained  in 
"teledb/emodels/reference"  is  called  "prem.ref.  dishask.csh  creates  a 
datafile  for  the  dishask  fortran  77  program  which  actually  does  the 
computations.  The  synthetic  Green’s  functions  will  be  put  in  the  directory 
"event/greens"  (e.g.,  yunnan/greens). 

Greens's  Function  Windowing  and  Decimation 

The  Green’s  functions  must  then  be  windowed  using  nearly  the  saine 
procedure  as  the  data.  Use  SAC  to  pick  the  onset  times  of  the  *.vss  Green’s 
functions  as  the  TO  marker.  (Only  use  the  *.vss  files.  The  other  *.vds  and 
*.clv)  files  will  be  automatically  windowed  by  SAC.)  Then  process  the  data 
using  the  window_greens.csh  script.  The  arguments  are  almost  identical  to 
that  in  window_data.csh.  The  windowed  Green’s  function  files  do  not  have 
any  special  suffix  appended  to  their  names,  but  are  simply  overwritten.  If  a 
mistake  is  made  here,  rerun  dishask.csh. 
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Moment  Tensor  Inversion 

At  last  the  Green's  functions  (for  one  source  depth)  and  the  data  have 
been  processed.  Inversion  can  now  proceed.  momvert__setiip.csh^  is 
provided  to  create  a  parameter  control  file  for  momvert,  the  inversion 
program.  This  is  another  somewhat  tedious  routine  which  asks  questions 
about  what  you  want  to  invert  for  and  so  on.  Finally,  momvert.csh  puts  it 
all  together  and  runs  the  momvert  program.  To  invert  for  other  source 
depths,  you  have  to  rerun  dishask.csh,  pick  TO  on  the  *.vss  Greens 
functions,  rerun  window  greens.csh  and,  possibly,  momvert_setup.csh. 

In  case  you  would  like  to  do  a  little  trial-and-error  modeling  of  the 
data,  a  C  Shell  script  called  synthetic.csh  has  been  written  to  run  the 
dishask  program  for  making  synthetics  for  particular  source  models.  I  r^ 
out  of  command  line  parameters  for  this  script  so  you  must  go  into  the  script 
and  edit  particular  entries  that  are  written  to  the  input  file  to  dishask.  Also, 
the  event.stachan  file  should  be  rebuilt,  preferably  with  another  nanae  to 
reflect  the  fact  that  synthetics  are  being  computed  and  not  Green's  functions. 
(The  mnvert  parameter  in  the  last  question  of  create_stachan.csh  must  be 
set  to  0,  not  1).  A  fast  way  to  make  an  appropriate  event.stachan  file  is  to 
edit  this  entry  in  a  text  editor. 

Random  Technical  Points 

momvert; 

The  arrays  in  momvert  are  dimensioned  for  a  total  of  4(XK)  data  points. 
For  example,  10  seismograms  which  are  each  400  points  long  could  be  used 
in  an  inversion.  Reducing  the  number  of  data  points  is  the  purpose  of 
waveform  decimation  and  the  option  in  momvert_setup.csh  for  choosing  an 
effective  decimation  in  the  inversion  program. 

20  time  function  elements  may  be  specified.  Time  function  element 
averaging  works;  use  an  odd  number  of  elements  if  you  want  to  average. 

Inversion  for  station  amplitude  factors  has  not  been  checked  and  is  not 
tmstworthy.  Do  not  use.  Assume  station  factors  of  unity. 

dishask: 

The  maximum  length  synthetic  time  series  allowed  is  2048  points. 
The  program  makes  sure  that  the  working  number  of  points  is  some  power 
of  2  for  FFT  use. 
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momplot:  .  u 

Source  code  for  this  program,  written  by  Jeff  Barker,  is  m  the 
teledb/src/momvert  directory.  It  is  an  interactive  program  which  computes 
synthetic  seismograms  from  the  Green's  functions  and  the  inversion  result 
and  plots  data  and  synthetic  around  moment  tensor  focal  mechanisms.  This 
program  will  be  available  once  I  repair  the  data  input  formats  and  find  the 
appropriate  graphics  libraries  on  our  new  system  at  Penn  State.  It  is  not 
available  for  the  short  course. 

JSPO  Softwsrc* 

The  DATASCOPE  software  is  available  from  the  IRIS/JSPC  via 
anonymous  ftp.  The  machine  address  is  jspc.colorado.edu.  You  will  need  a 
big  disk  to  hold  it. 
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Appendix 


This  appendix  contains  listings  of  the  C  Shell  scripts  which  drive 
Fortran  77  programs,  the  SAC  program,  and  the  evalresp  program, 
complete  listing  of  the  source  code  for  programs  written  here  is  available 


the  "teledb/src”  directory. 
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create_directory.csh 


pagel 


#!/bin/csh  -f 
# 

#  Create_directory.csh  -  shell  script  to  build  directories 

#  for  teleseismic  moment  tensor  inversion 

# 

if  ($#argv  !=  1)  then 

echo  "Usage:  $0  {directory  name}" 
exit 

endif 

# 

mkdir  ./$1 
mkdir  ./$1/data 
mkdir  ./$1/data/raw 
mkdir  ./$1 /data/process 
mkdir  ./$1 /data/invert 
mkdir  ./$1/sacplz 
mkdir  ./$1/synth 
mkdir  ./$1/greens 
mkdir  ./$1/invert 
mkdir  ./$1/db 

exit 
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create_stachan.csh 


pagel 


#!/bin/csh  -f 

create_stachan  .csh 

Set  up  input  file  for  the  sctable.f  program 
Run  sctable.f  to  produce  a  database  file  that  will 
be  used  in  the  dishask  program 

Note:  $PSUMOM  environmental  variable  is  needed  to 
find  all  directories,  sctable.f  will  also 
need  db.wfdisc  and  SAC  waveform  files 

if  ($#argv  !=  3  )  then 

echo  "Usage:  $0  {dbold  dbnew  eventdir}" 

echo  "Where," 

echo  "  dbold  =  working  database  table  name." 
echo  "  dbnew  =  new  database  table  name." 
echo  "  eventdir  =  event  directory  name." 
exit 
endif 

#  Set  up  directory  info  for  the  sctable  program 

set  dbpath  =  ($PSUMOM/$3/db) 

set  datapath  =  ($PSUMOM/$3/data) 

set  recpath  *  ($PSUMOM/emodels/receive) 

#  Write  file-pathnames  to  file  sctable.in 

echo  $dbpath  >  sctable.in 
echo  $datapath  »  sctable.in 
echo  $recpath  »  sctable.in 

echo  $1  »  sctable.in 
echo  $2  »  sctable.in 

#  Run  the  sctable  program  to  build  a  stachan  database  file  for 

#  each  waveform  of  interest 

sctable 
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dishask.csh 
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#!/bin/csh  -f 

dishask.csh 

Set  up  input  file  for  the  dishask.f  program 
Run  dishask.f  to  create  synthetic  Green's  functions 
for  use  in  the  inversion  program 

Note:  $PSUMOM  environmental  variable  is  needed  to 
find  ail  directories,  dishask.f  will  also 
need  db.wfdisc  and  db.stachan  files.  All 
Source,  receiver  and  reference  earth  model 
files  must  also  be  available  in  the  proper 
directories. 

if  ($#argv  !=  9  )  then 

echo  "Usage:  $0  {dbname  eventdir  refmodel  srcmodel  SACfile 
depth  tau  npts  dt}" 


echo  "  where  dbname 

=  working  database  name" 

echo  " 

eventdir 

=  event  directory  name." 

echo  " 

refmodel 

=  file  name  for  the  reference  earth 
model" 

echo  " 

srcmodel 

=  file  name  for  the  source  earth 
model" 

echo  " 

SACfil 

=  rootname  of  output  SAC  files" 

echo  " 

depth 

=  source  depth  (km)" 

echo  " 

tau 

=  time  width  of  triangular  time 
function  element" 

echo  " 

npts 

=  number  of  samples  in  the  synthetic 
time  series" 

echo  " 
echo"  " 

dt 

=  time  series  time  increment" 

echo  "  Note,  that  a  moment  of  1  .Oe+25  is  built  into  this  calculation." 
exit 

endif 

#  Write  SACfile  root  name  to  dishask.in 
echo  $5  >  dishask.in 

#  Set  up  directory  info  for  the  dishask  program 
set  dbpath  =  ($PSUMOM/$2/db) 

set  srcpath  =  ($PSUMOM/emodels/source/$4) 
set  recpath  =  ($PSUMOM/emodels/receive) 
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S0t  rofpath  =  {$PSUMOM/Grnod6ls/r8f/$3) 
set  grnpath  =  ($PSUMOM/$2/greens) 

#  Write  file-pathnames  to  file  dishask.in 

echo  $dbpath  »  dishask.in 
echo  $srcpath  »  dishask.in 
echo  $recpath  »  dishask.in 
echo  $refpath  »  dishask.in 
echo  $grnpath  »  dishask.in 

#  Write  database  name  to  file  dishask.in 
echo  $1  »  dishask.in 

#  Write  default  input  data  stream  to  dishask.in  for  computing 

’  #  Green's  functions 

echo  "1"  »  dishask.in 
echo  "  0.  0.  0."  »  dishask.in 

echo  "1.0  1.0  1.0  1.0  1.0  1 .0  "»  dishask.in 

echo  "1  .Oe+25"  $6  »  dishask.in 

echo  $7  "1"  »  dishask.in 

echo  "1 .0"  »  dishask.in 

echo  $8  $9  »  dishask.in 

#  Run  the  dishask  program  to  create  Green's  functions  for  this 
depth 

dishask  <  dishask.in 

#  All  Done! 
exit 
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#!/bin/csh  -f 

momvert.csh 

Run  the  Momvert  program. 

Note:  $PSUMOM  environmental  variable  is  needed  to 
find  all  directories,  gwindow.f  will  also 
need  the  db.stachan  file. 

if  ($#argv  !=  3  )  then 
echo  " " 

echo  "Usage:  $0  {dbname  eventdir  parmfile}" 
echo  "Where," 

echo  "  dbname  =  database  table  name  (as  in  dbname.waveparm)." 

echo  "  eventdir  =  event  directory  name." 

echo  "  parmfile  =  name  of  the  Momvert  parameter  file." 

exit 

endif 

set  invpath  =  ($PSUMOM/$2/invert) 
cd  $invpath 

#  Set  up  directory  info 
set  dbpath  *  ($PSUMOM/$2/db) 

#  Write  dbpath  to  file  momvert.in 
echo  $dbpath  >  momvert.in 

#  write  database  name  to  momvert.in 
echo  $1  »  momvert.in 

#  write  parameter  filename  to  momvert.in 
echo  $3  »  momvert.in 

#  Run  the  Momvert  program 
momvert  <  momvert.in 
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#!/bin/csh  -f 

momvert_setup.csh 

Set  up  parameter  files  and  waveform  info  files  for  the 
Momvert  program. 

Note:  $PSUMOM  environmental  variable  is  needed  to 
find  all  directories,  gwindow.f  will  also 
need  the  db.stachan  file. 

if  ($#argv  !=  4  )  then 
echo  -  " 

echo  "Usage:  $0  {dbname  eventdir  sacroot  parmfile} 
echo  "Where," 

echo  "  dbname  =  database  table  name  (as  in 

dbname.waveparm)." 

echo  "  eventdir  =  event  directory  name." 

echo  "  sacroot  =  Root  name  of  Green's  functions  Sac  files. 

echo  "  parmfile  =  name  of  the  Momvert  parameter  file." 

exit 

endif 

#  Set  up  directory  info 
set  dbpath  =  ($PSUMOM/$2/db) 

#  Write  Sacroot  to  file 
echo  $3  >  momparm.in 

#  Write  dbpath  to  file  momparm.in 
echo  $dbpath  »  momparm.in 

#  write  database  name  to  momparm.in 
echo  $1  »  momparm.in 

#  write  parmfile  to  file 
echo  $4  »  momparm.in 

#  Run  the  mvertw  program  to  create  a  dbname.waveparm  file 
mvertw  <  momparm.in 
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Run  the  momparm  program  to  create  a  Momvert  parameter  file 
Program  will  have  interactive  input. 

momparm 
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#!/bin/csh  -f 

prewave.csh 

Set  up  input  file  for  the  prewave.f  program 
Run  prewave.f  to  produce  new  database  files  and 
SAC  macro  file  for  instrument  response  correction 
in  addition  to  making  polezero  files  for  each 
channel 

Run  SAC  to  correct  for  instrument  responses 

Note:  $PSUMOM  environmental  variable  is  needed  to 
find  an  directories 

if  ($#argv  !=  5  )  then 

echo  "Usage:  $0  {dbold  dbnew  eventdir  refmodel  depth}" 
echo  "Where" 

echo  "  dbold  =  working  database  table  name." 

echo  "  dbnew  =  new  database  table  name." 

echo  "  eventdir  =  event  working  directory  name." 

echo  "  refmodel  =  spherical  earth  model  (without  .ref  suffix)." 

echo  "  depth  =  source  depth  for  travel  time  calculation."  exit 

endif 

#  Set  up  directory  info  for  the  prewave  program 

set  respath  =  ($PSUMOM/responses) 
set  plzpath  =  ($PSUMOM/$3/sacplz) 
set  dbpath  =  ($PSUMOM/$3/db) 
set  emodpath  =  ($PSUMOM/emodels) 
set  datapath  =  ($PSUMOM/$3/data) 

#  Write  to  file  prewave.in 

echo  $respath  >  prewave.in 
echo  $plzpath  »  prewave.in 
echo  $dbpath  »  prewave.in 
echo  $emodpath  »  prewave.in 
echo  $datapath  »  prewave.in 

#  write  database  names  to  prewave.in 

echo  $1  »  prewave.in 
echo  $2  »  prewave.in 
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#  specify  reference  earth  model 
echo  $4  »  prewave.in 

#  Specify  source  depth  for  travel  time  computation  here 
echo  $5  »  prewave.in 

#  Calculate  source  origin  time  and  put  into  prewave.in 
echo  $dbpath  >  rdorig.in 

echo  $1  »  rdorig.in 
set  torgin  =  'rdorigin  <  rdorig.in' 
set  torg  =  'epoch  $torgin' 
echo  $torg[2]  »  prewave.in 
echo  $torg[3]  »  prewave.in 
echo  $torg[4]  »  pr0wav0.in 

#  Run  the  prewave  program  to: 

#  Determine  source-receiver  distance 

#  Calculate-  P,  S  wave  travel  times 

#  Compute  polezero  file 

#  produce  new  wfdisc  file 

#  Produce  a  new  phase  associate  file 

#  Write  a  SAC  macro  for  instrument  correction 

prewave  <  prewave.in 

#  Run  the  plzero. script  in  the  sacpiz  directory 
cd  Splzpath 

chmod  755  plzero. script 
plzero. script 

#  Run  the  SAC  program  to  correct  data  for  instrument  response 

cd  $datapath/raw 
sac  <  sac.process.in 

#  Make  copies  of  the  old  database  files  into  new  database 

cp  $dbpath/$1  .origin  $dbpath/$2.origin 
cp  $dbpath/$1  .site  $dbpath/$2.site 

#  All  done! 
exit 
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#!/bin/csh  -f 
# 

#  rotate.csh 

#  Set  up  input  file  for  the  rotwave.f  program 

#  Run  rotwave.f  to  produce  new  database  files  and 

#  SAC  macro  file  for  vector  rotations 

#  Run  SAC  to  correct  for  instrument  responses 

# 

#  Note:  $PSUMOM  environmental  variable  is  needed  to 

#  find  all  directories,  rotwave.f  will  also 

#  need  db.wfdisc  and  db.origin  files  to  correct 

#  horizontal  data  channels. 

if  ($#argv  !=  3  )  then 

echo  "Usage:  $0  {dbold  dbnew  eventdir}" 
echo  "Where," 

echo  "  dbold  =  working  database  table  name." 
echo  "  dbnew  =  new  database  table  name." 
echo  "  eventdir  =  event  directory  name." 
exit 
endif 

#  Set  up  directory  info  for  the  rotwave  program 

set  dbpath  =  ($PSUMOM/$3/db) 
set  datapath  =  ($PSUMOM/$3/data) 

#  Write  to  file  rotwave.in 

echo  $dbpath  »  rotwave.in 
echo  $datapath  »  rotwave.in 

#  write  database  names  to  rotwave.in 

echo  $1  »  rotwave.in 
echo  $2  »  rotwave.in 

#  Run  the  rotwave  program  to: 

#  Determine  if  there  are  2  horizontal  components, per  station 

#  Check  SAC  header  for  sufficient  information  for  rotation 

#  Write  a  SAC  macro  for  data  rotation 

#  Produce  a  new  *.wfdisc  file  of  final  processed  data 
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rotwave  <  rotwave.in 

#  Run  the  SAC  program  to  correct  data,  if  needed,  and  then 

#  rotate 

cd  $datapath/process 
sac  <  sac.rotate.in 

#  Make  copies  of  the  old  database  files  into  new  database 

cp  $dbpath/$1  .origin  $dbpath/$2.origin 
cp  $dbpath/$1  .site  $dbpath/$2.site 

#  All  done! 
exit 
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#!/bin/csh  -f 


synthetic.csh 

Set  up  input  file  for  the  dishask.f  program 
Run  dishask.f  to  create  synthetic  seismograms. 

Edit  the  contents  of  this  script  to  input  appropriate 
source  parameters 

#  Note:  $PSUMOM  environmental  variable  is  needed  to 

#  find  all  directories,  dishask.f  will  also 

#  need  -db.wfdisc  and  db.stachan  files.  All 

#  Source,  receiver  and  reference  earth  model 

#  files  must  also  be  available  in  the  proper 

#  directories. 


if  ($#argv  !=  9  )  then 

echo  "Usage:  $0  {dbname  eventdir  refmodel  srcmodel  SACfile  depth 
tau  npts  dt}"  echo  "  where  dbname  =  working  database  name" 


echo 
echo 
echo 
echo 
echo 
echo 

echo 

echo 
echo 
echo 
echo  " 
synthetic, 
exit 


Note, 


eventdir  =  event  directory  name.’ 
refmodel  =  file  name  for  the  reference  earth  model" 
srcmodel  =  file  name  for  the  source  earth  model" 
SACfil  =  rootname  of  output  SAC  files" 
depth  =  source  depth  (km)" 
tau  =  time  width  of  triangular  time  function 
element" 

npts  =  number  of  samples  in  the  synthetic  time 
series" 

dt  =  time  series  time  increment" 

that  an  appropriate  dbname. stachan  file  must  exist  with' 
the  mnvert  parameter  set  to  0  to  calculate  a 


endif 

#  Write  SACfile  root  name  to  dishask.in 
echo  $5  >  dishask.in 

#  Set  up  directory  info  for  the  dishask  program 
set  dbpath  =  ($PSUMOM/$2/db) 

set  srcpath  =  ($PSUMOM/emodels/source/$4) 
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set  recpath  *  ($PSUMOM/emodels/receive) 
set  refpath  =  ($PSUMOM/emodels/ref/$3) 
set  grnpath  =  ($PSUMOM/$2/synth) 

#  Write  file-pathnames  to  file  dishask.in 

echo  $dbpath  »  dishask.in 

echo  $srcpath  »  dishask.in 

echo  $recpath  »  dishask.in 

echo  $refpath  »  dishask.in 

echo  $grnpath  »  dishask.in 

#  Write  database  name  to  file  dishask.in 

echo  $1  »  dishask.in 

#  Modify  these  input  data  streams  to  dishask.in  for  computing 

#  source  models  with  arbitrary  parameters. 

#  nfstyp  =0,  dislocation  parameters,  *1,  moment  tensor 

echo  "0"  »  dishask.in 

#  strike,  dip,  and  rake  in  degrees 

echo  "  30.  40.  60."  »  dishask.in 

#  moment  tensor  elements,  m11,m1 2, m1 3, m22,m23,m33 
echo  "1 .0  1.0  1.0  1 .0  1.0  1 .0  "  »  dishask.in 

#  seismic  moment  dyne/cm  and  depth 
echo  "1  .Oe+25"  $6  »  dishask.in 

#  tau  for  triangle  series  and  number  of  triangles  in  time  function 
echo  $7  "5"  »  dishask.in 

#  input  the  values  of  the  triangle  series,  the  area  of  the  time 

#  function  will  be  normalized  in  the  program, 
echo  "0.5  0.4  0.3  0.2  0.1"  »  dishask.in 

#  input  npts  and  delt-t  for  time  series 
echo  $8  $9  »  dishask.in 
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#  Run  the  dishask  program  to  calculate  synthetics 

dishask  <  dishask.in 

#  All  Done! 
exit 
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#!/bin/csh  -f 

window_data.csh 

Set  up  input  file  for  the  dwindow.f  program 
which  will  create  a  SAC  macro  file  for  windowing, 
filtering  and  decimation  of  waveform  data. 

Run  SAC  to  perform  this  stuff. 

Note:  $PSUMOM  environmental  variable  is  needed  to 
find  all  directories,  dwindow.f  will  also 
need  db.wfdisc  and  db.stachan  files. 

if  ($#argv  !=  7  )  then 
echo  " " 

echo  "Usage:  $0  {dbname  eventdir  corner  factor  twindow  decbh  decih}" 
echo  "Where," 
echo  " 
echo  " 
echo  " 
echo  " 

echo  " 

echo  " 

echo  " 
echo  " 

echo  " 

echo  " 

echo  " 

exit 
endif 

#  Set  up  directory  info  for  the  dwindow  program 

set  dbpath  =  ($PSUMOM/$2/db) 
set  datapath  =  ($PSUMOM/$2/data) 


dbnamG  =  database  table  name  (as  in  dbname.wfdisc). 

eventdir  =  event  directory  name.  Will  construct  name" 

PSUMOM/datadir/data." 

corner  =  corner  frequency  for  high  pass  butterworth 
filter." 

factor  =  constant  multiplicative  factor  to  convert 
units." 

twindow  =  data  waveform  time  window  wanted.  Will 
assume" 

5  seconds  before  tO  pick  time  as  start  time." 

decbh  =  decimation  factor  that  SAC  will  use  to 
decimate" 

the  20sps  data.  Use  a  value  of  1  for  no 
decimation." 

decih  =  decimation  factor  that  SAC  will  use  to 
decimate" 

the  1sps  data.  Use  a  value  of  1  for  no 
decimation." 
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it  Write  to  file  dwindow.in 

echo  $dbpath  >  dwindow.in 
echo  $datapath  »  dwindow.in 
it  write  database  names  to  dwindow.in 

echo  $1  »  dwindow.in 

#  Write  SAC  parameters  to  file,  also 

echo  $3  $4  $5  $6  $7  »  dwindow.in 

#  Run  the  dwindow  program  to  create  a  SAC  macro  file 
dwindow  <  dwindow.in 

#  Run  the  SAC  program  to  window  the  data 
cd  $datapath/process 

sac  <  sac.dwindow.in 

#  All  done! 
exit 
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#!/bin/csh  -f 

winclow_greens.csh 

Set  up  input  file  for  the  gwindow.f  program 
which  will  create  a  SAC  macro  file  for  windowing, 
filtering  and  decimation  of  waveform  Green's  functions. 

Run  SAC  to  perform  this  stuff. 

Note:  $PSUMOM  environmental  variable  is  needed  to 
find  all  directories,  gwindow.f  will  also 
need  the  db.stachan  file. 

if  ($#argv  !=  8  )  then 
echo  " " 

echo  "Usage:  $0  {dbname  eventdir  SACname  corner  factor  twindow 

decbh  decih}" 

echo  "Where," 

echo  "  dbname  =  database  table  name  (as  in  dbname. wfdisc)." 

echo  "  eventdir  =  event  directory  name.  Will  construct  name" 

echo "  PSUMOM/eventdir/greens." 

echo  "  SACname  =  Root  of  Green's  function  filenames  used  in" 

echo  "  Dishask  computations." 

echo  "  corner  =  corner  frequency  for  high  pass  butterworth 

filter." 

echo  "  factor  =  constant  multiplicative  factor  to  convert 

units." 

echo  "  twindow  =  data  waveform  time  window  wanted.  Will 

assume" 

echo  "  5  seconds  before  tO  pick  time  as  start 

time." 

echo  "  decbh  =  decimation  factor  that  SAC  will  use  to 

decimate  the  BH" 

echo  "  data.  Use  a  value  of  1  for  no  decimation." 

echo  "  decih  =  decimation  factor  that  SAC  will  use  to 

decimate  the  LH" 

echo  "  data.  Use  a  value  of  1  for  no  decimation." 

exit 

endif 

#  Input  SAC  root  name 
echo  $3  >  gwindow.in 
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#  Set  up  directory  info  for  the  dwindow  program 
set  dbpath  =  ($PSUMOM/$2/db) 

set  grnpath  =  ($PSUMOM/$2/greens) 

#  Write  to  file  gwindow.in 

echo  $dbpath  »  gwindow.in 
echo  $grnpath  »  gwindow.in 

#  write  database  name  to  gwindow.in 
echo  $1  »  gwindow.in 

#  Write  SAC  parameters  to  file,  also 
echo  $4  $5  $6  $7  $8»  gwindow.in 

#  Run  the  dwindow  program  to  create  a  SAC  macro  file 
gwindow  <  gwindow.in 

#  Run  the  SAC  program  to  window  the  Green's  functions 
cd  $grnpath 

sac  <  sac.g window. in 

#  All  done! 
exit 
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